home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
fortran
/
extenpr.zip
/
EXTEND.FOR
< prev
next >
Wrap
Text File
|
1988-06-04
|
16KB
|
709 lines
C
C This is a collection of subroutines that perform extended precision
C floating point arithmetic.
C
C By Mark P. Esplin
C 16 Shawsheen Rd.
C Billerica, Ma 01821
C
C Whole programs do not have to be written using extended
C precision. Only the part that needs extra precision need use it,
C then convert back to regular REAL*8 numbers.
C The routines for addition and multiplication use the same algorithms
C that are taught in grade school but using base 10000 numbers.
C This is not the most efficient way to do it, but
C it makes it possible to use standard FORTRAN77 (almost).
C These routines should work with any computer and compiler that supports
C INTEGER*2 and INTEGER*4 arithmetic.
C Each number is stored in an INTEGER*2 array of length N.
C The first word, IA(1), is the exponent in base 10000 digits.
C The least significant digit is in IA(2) and the most significant
C digit is in IA(N). If IA(N) > 5000 the number is negative.
C When IA(N) = 5000 the number is assumed invalid.
C Base 10000 digits have worse roundoff characteristics than binary.
C Since roundoff causes the last base 10000 digit to be inaccurate,
C that really means the last 4 decimal digits are inaccurate. Be
C sure to use more digits than is necessary to allow for roundoff.
C N must be at least 6 or some routines won't work correctly.
C
C Subroutines with their FORTRAN equivalent in this package are:
C
C ADD(A,B,C,N) C=A+B
C NEG(A,N) A=-A
C MUL(A,B,C,N) C=A*B
C INV(A,B,N,W) B=1/A W is work area of size 3*N
C EXSQRT(A,B,N,W) B=SQRT(A) W is work area of size 4*N
C EXHALF(A,N) A=1/2 much faster than using INV
C PRTEXT(A,N) WRITE(*,*) A
C WRTEXT(IU,A,N) WRITE(IU,*) A
C CVI2EX(I2,A,N) A=I2 I2 is INTEGER*2
C CVR8EX(R,A,N) A=R R is REAL*8
C CVEXR8(A,N,R) R=A R is REAL*8
C COPY(A,B,N) B=A
C IEXCMP(A,B,N,W) Compares two positive numbers A and B returns 1, 0, -1
C CHGPRE(A,NA,B,NB) Changes the precision of A form NA to NB
C
C There are other internal routines, but they would probably be of limited
C value for most users.
C
C
SUBROUTINE INCEX(IA,N)
INTEGER*2 IA(*),N
C This routine increments the mantissa of extended number IA.
C This routine is used by other routines to perform rounding.
INTEGER*2 IT,ICAR,I
ICAR=1
DO 10 I=2,N
IT=IA(I)+ICAR
ICAR=IT/10000
IA(I)=IT-ICAR*10000
C In almost all cases this routine will end after a few passes.
IF(ICAR.EQ.0) GOTO 20
10 CONTINUE
C Shift back if necessary
C The only time that a shift is necessary is when IA(N) = 5000 and
C the rest of the mantissa is zero.
20 CONTINUE
IF(IA(N).EQ.5000) THEN
IA(1)=IA(1)+1
DO 30 I=2,N-1
IA(I)=IA(I+1)
30 CONTINUE
IA(N)=0
C Don't need to round again since if a shift was necessary
C the least significant digit is zero.
ENDIF
RETURN
END
C
C
C
C
C
FUNCTION ICKCON(IA,N)
C This routine is used to check if an extended number IA is converging
C to an integer value.
C This routine would probably not be very useful for the average user.
C The value returned is 1 if IA has converged. The list significant digit
C is not checked since it might be just round off.
C This function return a 1 in IA has converged to an integer otherwise
C it returns a 0.
INTEGER*2 IA(*)
INTEGER*2 N,I
INTEGER*2 ITEST
ICKCON=0
C The number must have already converged past the N-2 th digit before
C this routine is called.
ITEST=IA(N-2)
DO 10 I=3,N-3
IF(IA(I).NE.ITEST) RETURN
10 CONTINUE
ICKCON=1
RETURN
END
SUBROUTINE NEG(IA,N)
C This routine changes the sign of extended precision number IA.
INTEGER*2 IA(*),N
INTEGER*2 I,J
J=2
DO 10 I=2,N
IF(IA(J).NE.0) GOTO 20
J=J+1
10 CONTINUE
C Number was zero
RETURN
C The first non-zero digit has been reached
20 IA(J)=10000-IA(J)
J=J+1
DO 30 I=J,N
IA(I)=9999-IA(I)
30 CONTINUE
RETURN
END
C
C
C
C
SUBROUTINE PRTEXT(IA,N)
C This routine writes extended precision numbers to output.
C This routine could use some work to make a more readable output.
INTEGER*2 IA(*),N
INTEGER*2 I,INEG
INEG=0
WRITE(*,*) ' '
IF(IA(N).GT.5000) THEN
WRITE(*,*) 'Minus'
CALL NEG(IA,N)
INEG=1
ENDIF
WRITE(*,900) IA(1)
900 FORMAT(' Exponent ',I6)
WRITE(*,910) (IA(I),I=N,2,-1)
910 FORMAT(10I5)
C Restore IA to previous value.
IF(INEG.NE.0) CALL NEG(IA,N)
RETURN
END
C
C
C
C
SUBROUTINE WRTEXT(IUNIT,IA,N)
C This routine writes extended precision numbers to file IUNIT.
C This routine could use some work to make a more readable output.
INTEGER*2 IUNIT,IA(*),N
INTEGER*2 I,INEG
INEG=0
WRITE(IUNIT,*) ' '
IF(IA(N).GT.5000) THEN
WRITE(IUNIT,*) 'Minus'
CALL NEG(IA,N)
INEG=1
ENDIF
WRITE(IUNIT,900) IA(1)
900 FORMAT(' Exponent ',I6)
WRITE(IUNIT,910) (IA(I),I=N,2,-1)
910 FORMAT(10I5)
C Restore IA to previous value.
IF(INEG.NE.0) CALL NEG(IA,N)
RETURN
END
C
C
C
C
C
SUBROUTINE NORMIZ(IA,N)
C This routine is used by other routines to normalize
C floating point number after addition and other such processes.
INTEGER*2 IA(*),N
INTEGER*2 I,ISHFT
IF(IA(N).GT.5000) THEN
C This is a negative number.
DO 10 I=N,2,-1
IF(IA(I).NE.9999) GOTO 20
10 CONTINUE
I=I+1
20 ISHFT=N-I
IF(IA(I).LE.5000) ISHFT=ISHFT-1
ELSE
C This is a positive number.
DO 30 I=N,2,-1
IF(IA(I).NE.0) GOTO 40
30 CONTINUE
C This number is zero
RETURN
40 ISHFT=N-I
IF(IA(I).GE.5000) ISHFT=ISHFT-1
ENDIF
IF(ISHFT.EQ.0) RETURN
IA(1)=IA(1)-ISHFT
C Shift left the appropriate amount.
DO 50 I=N,ISHFT+2,-1
IA(I)=IA(I-ISHFT)
50 CONTINUE
C Zero unknown digits.
DO 60 I=2,ISHFT+1
IA(I)=0
60 CONTINUE
RETURN
END
C
C
C
C
C
SUBROUTINE ADD(IA,IB,IC,N)
C This routine performs extended precision addition.
C IC = IA + IB. IA and IB are not changed.
INTEGER*2 IA(*),IB(*),IC(*)
INTEGER*2 N,IPA,IPB,NML,I
INTEGER*2 ISGNA,ISGNB,ISGNC,ISGEXT
INTEGER*2 IT,ICAR
ISGNA=1
IF(IA(N).GT.5000) ISGNA=-1
ISGNB=1
IF(IB(N).GT.5000) ISGNB=-1
C Figure out pointers for main loop
IF(IA(1).GE.IB(1)) THEN
IC(1)=IA(1)
IPA=2
IPB=IA(1)-IB(1)+2
NML=N-IPB+2
IF(NML.LT.2) THEN
C The difference between the two numbers is too large.
CALL COPY(IA,IC,N)
RETURN
ENDIF
IF(IPB.EQ.2) THEN
ICAR=0
ELSE
ICAR=(IB(IPB-1)+5000)/10000
ENDIF
ELSE
IC(1)=IB(1)
IPB=2
IPA=IB(1)-IA(1)+2
NML=N-IPA+2
IF(NML.LT.2) THEN
C The difference between the two numbers is too large.
CALL COPY(IB,IC,N)
RETURN
ENDIF
ICAR=(IA(IPA-1)+5000)/10000
ENDIF
C Main addition loop.
DO 10 I=2,NML
IT=IA(IPA)
IT=IT+IB(IPB)+ICAR
ICAR=IT/10000
IC(I)=IT-ICAR*10000
IPA=IPA+1
IPB=IPB+1
10 CONTINUE
C Do sign extended part.
IF(IA(1).GT.IB(1)) THEN
IF(ISGNB.GE.1) THEN
ISGEXT=0
ELSE
ISGEXT=9999
ENDIF
DO 20 I=NML+1,N
IT=IA(I)
IT=IT+ISGEXT+ICAR
ICAR=IT/10000
IC(I)=IT-ICAR*10000
20 CONTINUE
ELSE
IF(ISGNA.GE.1) THEN
ISGEXT=0
ELSE
ISGEXT=9999
ENDIF
DO 30 I=NML+1,N
IT=IB(I)
IT=IT+ISGEXT+ICAR
ICAR=IT/10000
IC(I)=IT-ICAR*10000
30 CONTINUE
ENDIF
C Shift right one place if it has overflowed during addition.
IF(ISGNA*ISGNB.EQ.1) THEN
ISGNC=1
IF(IC(N).GE.5000) ISGNC=-1
IF(ISGNA*ISGNC.EQ.-1) THEN
C An overflow has occurred.
IC(1)=IC(1)+1
ICAR=IC(2)
DO 40 I=2,N-1
IC(I)=IC(I+1)
40 CONTINUE
IC(N)=ISGEXT
C Do rounding and return.
IF(ICAR.GE.5000) CALL INCEX(IC,N)
RETURN
ENDIF
ENDIF
C Normalize answer.
CALL NORMIZ(IC,N)
RETURN
END
C
C
C
C
SUBROUTINE MUL(IA,IB,IC,N)
C This routine multiplies two extended precision numbers.
C N must be at least 4 or this routine won't work.
C IC = IA*IB
INTEGER*2 IA(*),IB(*),IC(*),N
INTEGER*2 ISNA,ISNB,IXC1,IXC2,ISHFT
C The computation is carried out with two extra placesIXC1 and IXC2
C This makes it so precision is not lost if the